## Load packages
pacman::p_load(
	broom, classInt, dismo, estimatr, Formula, fs,
	# ggeffects,
	ggplot2, geosphere, haven, labelled, lfe, lwgeom,
	magrittr, mgcv, needs, pacman, parallel,
	# pathological,
	pipeR, raster, RColorBrewer, readr, rgeos, reporttools,
	sf, sp, spdep, speedglm, stargazer, stringr, tidyverse, units
	)
prioritize(dplyr)
prioritize(purrr)
prioritize(lfe)

## ------------------------------------------------------------
## Projection setting
prj = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
## ------------------------------------------------------------

## --------------------------------------------------
## Functions
## --------------------------------------------------
extract_vec = magrittr::extract
extract_lst = magrittr::extract2
scale2 <- function(x, na.rm = TRUE) (x - mean(x, na.rm = na.rm)) / sd(x, na.rm)
log001p <- function(x) {
	x_logged = log(x + 0.01)
	return(x_logged)
	}
## Taken from https://www.statology.org/mode-in-r/
find_mode <- function(x) {
	u <- unique(x)
	tab <- tabulate(match(x, u))
	u[tab == max(tab)]
}
## Geodesic distance
# pointDistance = compiler::cmpfun(pointDistance)
all_distance = function(x, from_data, to_point) {
	dist = pointDistance(from_data[x,], to_point, type='GreatCircle', lonlat=TRUE)
	return(dist)
}
shortestDistance = function(x) {
	dist = min( pointDistance(fromDat[x,], toDat, type='GreatCircle', lonlat=TRUE) )
	return(dist)
}
shortestDistance_rev = function(x, from_data, to_point) {
	dist = min( pointDistance(from_data[x,], to_point, type='GreatCircle', lonlat=TRUE) )
	return(dist)
}
meanDistance = function(x) {
	dist = mean( pointDistance(fromDat[x,], toDat, type='GreatCircle', lonlat=TRUE) )
	return(dist)
}
## --------------------------------------------------
## Vector/list of variable names to formula obj
chr2fml <- function(dv, idv_list) {
	idv = sapply(idv_list, str_c, collapse = " + ")
	idv = idv[str_count(idv) != 0]
	idv = str_c(idv, collapse = " + ")
	mod_obj = as.formula( str_c(dv, "~", idv) )
	return(mod_obj)
}
## --------------------------------------------------
## Extract xvar from fomula object
fm2xvar <- function(fm) {
	vars = strsplit(as.character(fm), split = " ~ ", fixed = TRUE)[[3]]
	vars = unlist(strsplit(vars, split=" + ", fixed = TRUE))
	vars = unlist(strsplit(vars, split=" * ", fixed = TRUE))
	vars = unique(vars)
	vars = str_replace_all(vars, "\n  ", "")
	return(vars)
}
## --------------------------------------------------
fm2xvarWspline <- function(fm) {
	vars = strsplit(gsub("as.factor\\(|s\\(|\\)\\(|te\\(|\\)","", fm, perl = TRUE), split=" + ",fixed=TRUE)[[3]]
	vars = unlist(strsplit(vars, ","))
	vars = str_trim(vars, side="both")
	return(vars)
}
## --------------------------------------------------
metallic_ratio <- function(baseline, n = 1) {
	out = (n+sqrt(n^2+4))/2
	out = baseline*out
	return(out)
}
gld_ratio <- function(baseline) {
	out = baseline*(1+sqrt(5))/2
	return(out)
}
slv_ratio <- function(baseline) {
	out = baseline*(1+sqrt(2))
	return(out)
}
plt_ratio <- function(baseline) {
	out = baseline*sqrt(3)
	return(out)
}
## --------------------------------------------------
scalebar_rev <- function(loc, length, unit="km", division.cex=.7, line_width = 12,...) {
	if(missing(loc)) stop("loc is missing")
	if(missing(length)) stop("length is missing")
	x <- c(0,length/c(4,2,4/3,1),length*1.1)+loc[1]
	# y <- c(0,length/(line_width*3:1))+loc[2]
	y <- c(0,length/(10*3:1))+loc[2]
	cols <- rep(c("black","white"),2)
	for (i in 1:4) rect(x[i],y[1],x[i+1],y[2],col=cols[i])
	for (i in 1:5) segments(x[i],y[2],x[i],y[3])
	labels <- x[c(1,3)]-loc[1]
	labels <- append(round(d2km(labels),0),paste(round(d2km(x[5]-loc[1]),0),unit))
	text(x[c(1,3,5)],y[4],labels=labels,adj=.5, cex=division.cex, offset=.1,pos=3)
}
d2km <- function(d){
	out <- d*60*1.852
	return(out)
}
km2d <- function(km){
	out <- (km/1.852)/60
	return(out)
}
## --------------------------------------------------
## Rescale to [0,1] range
rescaleIt01 <- function(x, xrange = range(x, na.rm = TRUE)) {
	x2 = (x - xrange[1])/(xrange[2] - xrange[1])    ## range [0,1]
	x2 = x2
	return(x2)
}
## Centering function
centering <- function(x, na_rm = TRUE) {
	x = x - mean(x, na.rm = na_rm)
	return(x)
}
## --------------------------------------------------
## Generate convex hull from raster input
raster2convex_hull <- function(raster_obj, prj) {
	## ------------------------------------------------------------
	## Extract raster boundaries
	tmp_boundry = raster_obj
	tmp_boundry[tmp_boundry == 0] = NA
	tmp_boundry = boundaries(tmp_boundry, asNA = TRUE, directions = 8)
	## ------------------------------------------------------------
	## Cropping window polygon
	tmp_grd = as(tmp_boundry, "SpatialGridDataFrame")
	names(tmp_grd) = "raster_val"
	## ------------------------------------------------------------
	## Extract grid coordinates
	cods = coordinates(tmp_grd)  %>%
		tbl_df  %>%
		mutate(
			raster_val = as.vector(tmp_grd@data[,1])
			)   %>%
		filter(!is.na(raster_val))  %>%
		rename(x_cord = s1, y_cord = s2)
	## ------------------------------------------------------------
	## Create shapefile window: Convex hull
	cods = cbind(cods$x_cord,cods$y_cord)
	polyz = convHull(cods)
	polyz = polygons(polyz)
	st_pol = st_as_sf(polyz)
	st_crs(st_pol) = prj
	## ------------------------------------------------------------
	return(st_pol)
}
## --------------------------------------------------
# generate_polynomial <- function(x, degree, var_name = NULL){
#     output_matrix = matrix(nrow = length(x), ncol = degree)
#     for (i in 1:degree) {
#         output_matrix[,i] = x^i
#     }
#     if (!is.null(var_name)) colnames(output_matrix) = str_c(var_name, c("", 2:degree))
#     return(output_matrix)
# }
polynomial_varname <- function(var_vec, degree) {
	for (i in 1:length(var_vec)) {
		tmp_poly = str_c( "I(", var_vec[i], "^", 2:degree, ")")
		if (i > 1)  {
			output = c(output, tmp_poly)
		}   else    {
			output = tmp_poly
		}
	}
	output = c(var_vec, output)
	return(output)
}
## ----------------------------------------------------------------------
## Vector/list of variable names to formula obj
chr2fml_felm <- function(dv, idv_list, fe_list = NULL, endog_var = NULL, instrument_var = NULL, se_cluster_list = NULL) {
	require(Formula)
	## IDVs
	idv = sapply(idv_list, str_c, collapse = " + ")
	idv = idv[str_count(idv) != 0]
	idv = str_c(idv, collapse = " + ")
	## FEs
	if (is.null(fe_list)) {
		fixed_effects = "0"
	}	else	{
		fixed_effects = sapply(fe_list, str_c, collapse = " + ")
		fixed_effects = fixed_effects[str_count(fixed_effects) != 0]
		fixed_effects = str_c(fixed_effects, collapse = " + ")
	}

	## Cluster
	if (is.null(se_cluster_list)) {
		se_cluster = "0"
	}	else	{
		se_cluster = sapply(se_cluster_list, str_c, collapse = " + ")
		se_cluster = se_cluster[str_count(se_cluster) != 0]
		se_cluster = str_c(se_cluster, collapse = " + ")
	}

	## Convert into
	mod_obj = as.Formula( str_c(dv, "~", idv, "|", fixed_effects, "|0|", se_cluster) )
	if (!is.null(instrument_var)) {
		mod_obj = as.Formula( str_c(dv, "~", idv, "|", fixed_effects, "| (", endog_var, "~", instrument_var, ") |", se_cluster) )
	}
	return(mod_obj)
}
## ----------------------------------------------------------------------
model_info_lm <- function(model, x2extract, use_this_se = NULL, alpha = 0.95) {
	## Model info
	tmp_summary = summary(model)
	tmp_coef = tmp_summary$coefficients[x2extract, 1]
	## SE
	if (is.null(use_this_se)) {
		tmp_se = tmp_summary$coefficients[x2extract, 2]
	}   else    {
		tmp_se = use_this_se
	}
	## DF
	model_df = tmp_summary$df[2]
	if (is.na(model_df)) model_df = tmp_summary$rdf ## speedlm
	## CI
	qpoint = 1 - (1-alpha)/2
	crit_val = qt(qpoint, df = model_df)
	ci_lower = tmp_coef - crit_val*tmp_se  ## lower CI
	ci_upper = tmp_coef + crit_val*tmp_se  ## upper CI
	## Output vector
	output = c(tmp_coef, tmp_se, ci_lower, ci_upper)
	names(output) = c("coef", "se", "ci_lower", "ci_upper")
	return(output)
}
## --------------------------------------------------
model_info_gam <- function(model, x2extract, alpha = 0.95) {
	x2extract = str_replace_all(x2extract, "\\(", "\\\\(")
	x2extract = str_replace_all(x2extract, "\\)", "\\\\)")
	## Model info
	tmp_summary = summary(model)
	tmp_coef = tmp_summary[1]$p.coeff[str_detect(names(tmp_summary[1]$p.coeff), x2extract)]
	tmp_se = tmp_summary[2]$se[str_detect(names(tmp_summary[2]$se), x2extract)]
	model_df = tmp_summary$residual.df
	## CI
	qpoint = 1 - (1-alpha)/2
	crit_val = qt(qpoint, df = model_df)
	ci_lower = tmp_coef - crit_val*tmp_se  ## lower CI
	ci_upper = tmp_coef + crit_val*tmp_se  ## upper CI
	## Output vector
	if (length(tmp_coef) == 1) {
		output = c(tmp_coef, tmp_se, ci_lower, ci_upper)
		names(output) = c("coef", "se", "ci_lower", "ci_upper")
	}   else    {
		output = cbind(tmp_coef, tmp_se, ci_lower, ci_upper)
		colnames(output) = c("coef", "se", "ci_lower", "ci_upper")
	}
	return(output)
}
## --------------------------------------------------
model_info_robustlm <- function(model, x2extract, alpha = 0.95) {
	require(broom)
	x2extract = str_replace_all(x2extract, "\\(", "\\\\(")
	x2extract = str_replace_all(x2extract, "\\)", "\\\\)")
	tidied = tidy(model, conf.level = alpha)    %>%
		filter(str_detect(term, x2extract)) %>%
		# filter(term %in% x2extract) %>%
		select(
			coef = estimate, se = std.error,
			ci_lower = conf.low, ci_upper = conf.high
			)   %>%
		as.matrix
	return(tidied)
}
## --------------------------------------------------
predict_GAM <- function(mod, var, vals, quadratic=FALSE, fe_list=NULL){

	## Prepare design matrix
	vars <- names(mod$model)
	n <- length(vals)
	X <- as.data.frame(matrix(NA,nrow=n,ncol=length(vars)))
	names(X) <- vars
	vars. <- vars[!(grepl("as.factor",vars,fixed=T)|grepl("Intercept",vars,fixed=T))]
	X[,vars.] <- as.data.frame(matrix(apply(mod$model[,vars.],2,function(x){median(x,na.rm=T)}),nrow=n,ncol=length(vars.),byrow=T))
	X[,var] <- vals     ## main predictor
	for(j in 1:length(vars.)){X[,vars.[j]] <- as.numeric(as.character(X[,vars.[j]]))}

	## Quadratic terms
	if (quadratic) {
		quadz <- names(mod$model)[grep("I(",names(mod$model),fixed=TRUE)]
		for(k in 1:length(quadz)){
			quadd <- names(mod$model)[apply(as.data.frame(names(mod$model)),1,function(x){grepl(x,quadz[k])})]
			X[,quadz[k]] <- X[,quadd]^2
		}
	}

	## Fixed effects
	if (length(fe_list)>0) {
		for(k in 1:length(fe_list)){
			X[,which(grepl(names(fe_list[k]),vars))] <- as.factor(fe_list[[k]])
		}
	}

	## Draw predicted values
	X = as.data.frame(X)
	names(X) = gsub(")","",gsub("as.factor(","",names(X),fixed=TRUE))
	preds = predict.gam(mod, newdata=X, type="link", se.fit=TRUE)
	preds = cbind(preds$fit, preds$fit-1.96*preds$se.fit, preds$fit+1.96*preds$se.fit)
	colnames(preds) = c("y_hat", "ci_lower", "ci_upper")
	preds = mod$family$linkinv(preds)  ## (inverse) link function

	return(preds)
}


## EOF
